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1 Introduction: “converting chemical information into a ge¬ 
ometrical form” 

Alan Turing’s celebrated paper “The Chemical Basis of Morphogenesis” |22] pub¬ 
lished in 1952 in the Philosophical Transactions of the Royal Society of London is 
a typical example of a pioneering and inspired work in the domain of mathematical 
modelling. 

1. The paper presents a new key idea for solving an old problem. As Lionel 
Harrison said in his 1987 paper [0] “it is a theoretical preconception proceeding 
experience”. 

2. It contains right away the germ of quite the whole theory associated to the 
new ideas. 

3. Its forsightedness is striking. It anticipates by many years its experimental 
conhrmations and mathematical developments. 

The key idea is formulated from the outset in the first sentence: 

“It is suggested that a system of chimical substances, called morphogens, 
reacting together and diffusing through a tissue, is adequate to account 
for the main phenomena of morphogenesis.” 
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Later (1953), Turing used a striking formulation (see below figure [^: 

“It was suggested in Turing (1952) that this might be the main means 
by which the chemical information contained in the genes was converted 
into a geometrical form.” 

Turing’s colleague Claude Wilson Wardlaw, a botanist at the Department of 
Cryptogamic Botany at Manchester University who wrote with Turing the paper 
“A diffusion reaction theory of morphogenesis in plants” said in 1952 that 
Turing’s working hypothesis was that: 

“a localized accumulation of gene-determined substances may be an es¬ 
sential prior condition [for cell differentiation].” (p. 40) 

and that what is needed in addition to biochemistry for understanding emerging 
spatial forms in morphogenetic processes is a “patternized distribution of morpho¬ 
genetic substances”. There exist many homologies of organization between different 
biological species and it must therefore exist general morphogenetics mechanisms 
largely independent of specihc genes. As he claimed, “certain physical processes are 
of very general occurrence” (p. 46). 

Wardlaw summarizes Turing’s key idea in the following way: 

“In an embryonic tissue in which the metabolic substances may initially 
be distributed in a homogeneous manner, a regular, patternized distri¬ 
bution of specihc metabolites may eventually result, thus affording the 
basis for the inception of morphological or histological patterns.” (p. 44) 

So, for understanding how a spatial order can emerge from biochemical reactions 
genetically controlled - which, according to Turing, is the main problem of morpho¬ 
genesis - Turing dehnes from the outset a form as a breaking of homogeneity of some 
spatially extended biological tissue and, therefore, as a breaking of the symmetries 
underlying such an homogeneity. 

2 The framework 

2.1 Turing diffusion-driven instability 

The components of Turing’s theorization are the following. Inside the biological tis¬ 
sue under consideration, chimical reactions occur, which we will call internal chimical 
dynamics. In the spatial extension of this substrate, processes of spatial diffusion 
occur, which we will call external spatial dynamics. The key idea is that the coupling 
of these two very different kinds of dynamics can trigger, under certain conditions, 
morphogenetic processes, which can be mathematically modelled by using what are 
called since Turing reaction-diffusion differential eguations. Why and how? Because 
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the external spatial diffusion can destabilize, under certain conditions, the internal 
chimical equilibria. 

We must emphasize the fact that the notion of diffusion-driven instabilities is to 
some extent paradoxical. Indeed, diffusion is a stabilizing process and therefore the 
idea amounts to posit that the coupling of two stabilities can induce an instability! 

We will see that Turing assumes that there exists only one equilibrium of the 
internal chimical dynamics. A diffusion-induced instability could therefore make the 
chimical state diverge, but in general non-linearities of the equations bound such 
divergences. Another possibility would be that there exist several equilibria. Then 
Turing instabilities would induce bifurcations from one equilibrium state to another 
one. It is this idea that has been worked out in the late 1960s by Rene Thom ra. 
[2T] to explain also morphogenesis. 

2.2 Turing’s objective 

Today, reaction-diffusion equations are mainly used to explain the formation of pat¬ 
terns in material substrates. But Turing’s objective was deeper and more ambitious 
and concerned embryogenesis. As he claimed in his Introduction, 

“The purpose of this paper is to discuss a possible mechanism by which 
the genes of a zygote may determine the anatomical structure of the 
resulting organism.” (p. 37) 

As this general objective was too ambitious, he assumed many simplifications. 
The first simplification was to eliminate any direct reference to specific genes and to 
reduce the internal chemical dynamics to reaction equations between concentrations 
of morphogens. As explains Philip Maini in [9] , a morphogen is 

“a chemical to which cells respond by differentiating in a concentration- 
dependent way.” 

Turing was inspired by what Waddington [23] called “form producers” or “evoca¬ 
tors”. Morphogens are controlled by genes which catalize their production, but, 
contrary to genes, they can diffuse in the developing tissues and carry the positional 
information (e.g. in the sense of Wolpert) which is needed for morphogenesis]^ 

The second simplification made by Turing was to eliminate the mechano-chemical 
aspects of embryogenesis, although he was aware of their importance during the 
development. These aspects will be worked out later by specialists such as George 
Oster and James Murray 

But, even so drastically simplified, the search for good mathematical models 
remains a 

“a problem of formidable mathematical complexity” (p. 38). 

^For an introduction to the concept of “positional information” in Waddington, Wolpert, Good¬ 
win and Thom, see Petitot |18j . 
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Figure 1; Turing’s general equation for morphogenesis. 
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Figure 2: Turing’s anticipation of the richness of the general equation for morpho¬ 
genesis. 


In fact, it seems that Turing was looking for a kind of universal equation for 
morphogenesis. In his last paper “Morphogen theory of phyllotaxis”, which remained 
unpublished because of his suicide and is kept at the King’s College Archives, he 
proposed the equation 


+ /m (bl, • • ■ , Fm) 

where the F^ are the respective concentrations of the M morphogens, the spatial 
Laplacian, that is the diffusion operator, the /i^ the diffusibility coefficients, and 
the fm the reaction equations (see figure 

The point was that, when you vary the fm and the /i^, the solutions of such a 
universal equation can be extremely diverse (see figure]^. Hence the idea that, as 
with Newton’s equation for Mechanics, it could be possible to classify a lot of very 
different kinds of forms using the same general equation. 
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Indeed, we can vary three classes of parameters: 

1. For the chemical part, the eigenvalues provided by the spectral analysis of the 
linearized system of the fm- 

2. The diffusibility coefficients /i^. 

3. For the geometrical part, the eigenfunctions of the Laplacian operator (har¬ 
monic analysis). 

2.3 Turing’s foresightedness 

In his 1990 paper “Turing’s theory of morphogenesis. Its influence on modelling 
biological pattern and form” na, James Murray claims that Turing’s 1952 paper is 
“one of the most important papers in theoretical biology of this century” (p. 119). 
Indeed, 

“What is astonishing about Turing’s seminal paper is that, with very 
few exceptions, it took the mathematical world more than 20 years to 
realise the wealth of fascinating problems posed by his theory. What is 
even more astonishing is that it was closer to 30 years before a signihcant 
number of experimental biologists took serious notice of its implications 
and potential applications in developmental biology, ecology and epi¬ 
demiology.” (p. 121) 

These inspired anticipations proved to be exact in chemistry. There exist today a 
lot of models of chimical reaction-diffusion phenomena: clocks, travelling waves, etc. 
Their analysis constitutes a rapidly expanding research domain while, at Turing’s 
time, no empirical example was known. Turing discovered theoretically the basic 
phenomenon and was the hrst to compute simulations on the computer he had 
himself constructed at Manchester. In embryogenesis, the exact limits of validity of 
Turing model are still under discussion. 

3 The context 

3.1 The bibliography 

It is interesting to look at Turing’s bibliography, which is very short. First, it in¬ 
cludes two books which are not really used, the Theory of Elasticity and Magnetism 
of James Jeans (1927) and The permeability of natural membranes of Hugh Dawson 
and James Danielli (1943). Then, it cites a fundamental paper of Leonor Michaelis 
and Maud Menten (1913) on Die Kinetik der Invertinwirkung [T3] whose pioneering 
mathematical model is typical of the internal chimical dynamics used by Turing. 
Finally, there are three masterpieces on embryology and morphogenesis: Charles 
Manning Child’s “summa” Patterns and problems of development (1941), SirD’Arcy 
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Thompson’s masterpiece On Growth and Form (1942) [T] and Conrad Hal Wadding- 
ton’s key work on Organizers and Genes (1940) [2S]. Introduced by Hans Spemann, 
“organizers” were thought to be the cause of the embryological induction observed 
when the tissues of some part of an embryo (e.g. a leg) were transplanted in another 
part (e.g. the head). The idea (much speculative at that time) was that there must 
exist chemical signals triggering cellular differenciations. It is in the second part of 
this work that Waddington assumed that, through morphogens, gene concentrations 
could be important for cellular differenciation and that the developmental units of 
an organism are “morphogenetic helds”. 

3.2 The kinetic model 

Building on previous very precise numerical experimental data collected by Victor 
Henri (1903), the Michaelis-Menten model of the kinetics of invertase enzyme (1913) 
was the hrst to explain the catalysis of the hydrolysis of sucrose into glucose and 
fructose. Let E be an enzyme bounding with a substrate S to give a complex ES 
which converts itself into a product P through a chain of two elementary chemical 
reactions: 


E + s'^ES^ E + P 

k2 

where the ki are the rate constants of the reactions. Let us denote by [X] the 
concentration of X. Then the law of mass action saying that a reaction rate is 
proportional to the product of the concentrations of the reactants implies the system 
of nonlinear differential equations]^ 


[^] = -fci [E] [^] + k2 
^ [E] = -k, [E] [^] + k2 [i?^] + k, 

= fci [E] [^] - k2 - k3 

. IF] = h [BS] 

where the relation [E] -|- [ES] = 0 implies the conservation law [E] + = Eq = 

constant. Under an hypothesis of adiabaticity according to which the equilibrium 

between S and ES is “instantaneous”, that is [S'] = 0, then ki [E] [S'] = ^2 [ES], 
|BS| = I [E] [SJ = I [S] (B„ - [ES]), [BSI (l + | [S]) = | [SJ E„ and 


[ES] = I [S] 





{^ + ^ [“S']) 


with K = and therefore 

ki ’ 

is the traditional notation for the temporal derivative 


Eo[S] 

K+[S] 
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[^] = h 


Eo[S] 

K+[S] 


Turing’s numerical example 


So, Turing start with morphogens diffusing and reacting inside a tissue. Diffusion 
flows from regions of strong concentrations towards regions of weak concentrations 
with a velocity proportional to the gradients of the concentrations and to the dif- 
fusibility coefficients. According to the law of mass action, reaction rates are propor¬ 
tional to the product of concentrations. Hence a huge variety of nonlinear differential 
equations. 

Turing gives several examples and develop one of them in minute detail in his 
§10 “A numerical example”. He considers a ring of A^ = 20 cells and two mor¬ 
phogens X and Y and makes several numerical assumptions on their size, the dif- 
fusibility constants, the permeability of membranes (it is here that the reference to 
Dawson-Danielli is used), etc. With great acuity, he argues that the system must 
be (thermodynamically) open and include a “fuel substance” A providing it with 
energy through its degradation into another substance B. 


“In order to maintain the wave pattern a continual supply of free energy 
is required. It is clear that this must be so since there is a continual 
degradation of energy through diffusion. This energy is snpplied throngh 
the ‘fnel snbstances’ (A, B in the last example), which are degraded into 
‘waste products’.” (p. 65) 

For modelling catalysis, Turing add three other snbstances C, C", and W. 

It must be emphasized that it will be only thirty years later that open thermo¬ 
dynamical systems ont of equilibrinm will be systematically investigated (see e.g. 

Prigogine’s dissipative structures). 

4.1 The internal chimical dynamics 

The system of 7 elementary reactions proposed by Turing is the following: 


[1] F + X ^ IF 

[2] W + A ^2Y + B 


iXF 


34 = T 
^ 16 


rate: 

(instantly) 

A = 1000 = constant (fnel snbstance) 
rate: 
rate: 
rate: 

(instantly) 

c = C' = 10-3 (1 + 7) 

rate: §1036"= §(1 + 7 ) 

So, X converts into Y at the rate A [50XF + — 55 (1 + 7 )] (because of [1], 

[3], and [7]) while self-reprodncing (because of [4]) at the constant rate A 


[3] 2X ^ IF 

[4] A ^ X 
[b]Y ^B 

[Q]Y + C 

[7]C' ^X + C 


7 

64 A 

fIO- 
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destroying V (because of [5]) at the rate ^Y. So the kinetic equations for the 
time varying concentrations X [t), Y [t) of the morphogens X,Y- i.e. the internal 
chimical dynamics - are 


\X = X [-50XY - 7X2 + 57 + 55 ^] = f (x, Y) 

\ y = X [5oxy + 7X2 _ 55 _ 55 ^ _2Y]=g (X, Y) 

where 57 in the hrst equation is 55 + 2 with 2 coming from [4] and —2y in the 
second equation comes from [5]. 

4.2 Equilibria and linearization 

A chemical internal equilibrium corresponds to values (Xe, 17) such that / (Xe, 17) = 
0 and g (Xe, 17) = 0. For 7 = 0 the system is 


f X = X [-50Xy - 7X2 + 57] = / (X, Y) 

\ y = X [5oxy + 7X2 _ 55 _ 2 y] = ^ (X, y) 

and an evident equilibrium is Xq = Yq = 1. Another is Xi = — y,yi = 1 but it is 
non physical since a concentration cannot be negative. 

Then, Turing linearizes the system near the equilibrium (Xq, Tq) and analyzes 
the stability of the linear system. The method was already well known at his time: 
the matrix of the linearized system is the Jacobian Jq of {f,g} at (Xq, To) and the 
stability depends upon the fact that the real part of all eigenvalues A of Jq are < 0. 
Let us briefly remind non mathematicians of it. It is immediate to compute Jq for 
7 = 0 : 


and, as 


we get 


Jn = 


A/ dl 

■F ¥ 

og 09 
ax av 


f § = « = -I [-50>o - 14Xo] = -2 
)w = b=Y2 [- 50 ^ 0 ] = -i = -1.5625 
1 S = c = M [50>o X 14Xo] = 2 
l§ = ^=M[50^o-2] = i = i = 1.5 



In a small neighbourhood of (Xq, Yq) we can write X = Xq + x, Y = YQ + y and 
write at hrst order 


X 

y 



i.e. 


X = ax + by 
y = cx + dy 



As the system is linear, we look at solutions of the form 


_ pAt f xo 

y{t)) V^o 

where (xo,yo) is the state of the system at time t = 0. They are straight trajecto¬ 
ries on the line (0,0) — (xo,yo) with an exponential temporal law. Computing the 
derivatives in two different ways, we get 



that is an equation linking A to Jo: 


Ae^* 



A 



(^0 - V) = 0 

If (a;o,|/o) 7^ (0,0), then (x,y) ^ (0,0) and this equation can be satished only if the 
determinant Det (Jq — AJ) vanishes. The equation Det (Jq — AJ) = 0 is called the 
eharaeteristie equation of the linear system. It is a polynomial equation of degree 2 
which writes 


A^ 

A^ 


Det 


(a-\ b \ 
V c d-\) 


{a + d) \ + ad — be 
Tr ( Jq) a -|- Det (Jo) 
X^-SX + P 


0 

0 

0 

0 


where the sum of diagonal terms Tr ( Jq) = a -|- J = S', called the trace of the 
matrix Jq, gives the sum S of the solutions and the determinant Det (Jo) of Jq 
gives their product P. As the discriminant of the equation is A = S*^ — 4P = 
Tr (Jo)^ — 4 Det (Jo), the solutions are given by the well known formula 


a.A(s±7a) 

A± = i (^Tr (Jo) ± \/Tr(Jo)"-4Det(J„)^ 

and any solution of the linear system is a linear combination of the two solutions 
with A±. 

Let us suppose now that X = a + iu has real part Re (A) = a and imaginary 
part Im (A) = u. Then, is an oscillation modulated by 

the real exponential e“*. If a > 0, diverges exponentially when t —)■ -|-cx) and 
the corresponding trajectories go to inhnity and are unstable. On the contrary, if 
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a < 0 , e"* converges exponentially towards 0 when t —)■ +cx) and the corresponding 
trajectories go to equilibrium and are stable. In Turing’s example, 


Tr (Jo) — a + d — —2 2 ~ ~2 

3 / 25\ 1 

Det {Jo) = ad - 6 c = -2 x - - ( j ^ ^ = 3 

A± = --(l±i), Re(A±) = --<0 


and the two eigenvalues have negative real parts. The internal chemical equilibrium 
(Xo,To) (i-e. (a^o, Jo) = 0) is therefore stable. 

More generally, in the 2-dimensional case, the system is stable if and only if 


f Tr (Jo) = a -|- J < 0, here — | < 0 
( Det (Jo) = ad — bo 0, here | > 0 

Indeed, we must have Re (A±) <0. If A = Tr (Jo)^ — 4 Det (Jo) < 0, then A+ and 
A_ are complex conjugate eigenvalues and we must have Tr (Jo) < 0, and of course 
Det (Jo) > 0 because otherwise we would have A > 0. If A > 0, then A+ and A_ are 
real eigenvalues and they must be both < 0. This imply that the greatest eigenvalue, 
namely A+ = Tr (Jo)-I-V^ must be < 0, which implies Tr (Jo) < — VA. So Tr (Jo) < 
0, and, as Tr (Jo)^ > A = Tr (Jo)^ — 4 Det (Jo), we have also Det (Jo) > 0. 

5 Diffusion-driven instability 

After having dehned the internal chemical equilibrium and analyzed its stability, 
Turing explains how a spatial diffusion of the morphogens X, Y can induce an 
instability. Let us summarize his computations. 

5.1 The reaction-diffusion model 

Let r = I,-- - ,N label the positions of the N cells in the ring. Concentrations 
X, Y are then functions X{r,t) and Y{r,t) of time t and spatial position r. In a 
continuous model, the spatial positions would be parametrized by an angle 6 ^ G 
and concentrations would be functions X{9, t) and Y{9, t) (angles 9r = 27r^ retrieve 
the discrete case). We start with an homogeneous initial state where X{r,t) = Xq 
and Y (r, t) = To everywhere and we apply diffusion using the Laplace operator 
A = and its discrete approximation XF (r) = F {r — 1) — 2F (r) -|- F (r -|- 1) for 
any function F (r). 

We get that way the reaetion-diffusion equations 


'A(r,t) = /(A(r,t),F(r,t)) + 

p(A(r - l,t) -2X{r,t) + X{r + l,t)) 

Y{r,t) = g{X{r,t),Y{r,t)) + 

^ iy{Y{r-l,t)-2Y{r,t)+ Y{r+ l,t)) 
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where /i and v are the respective coefficients of diffusibility of X and Y. Some 
technical aspects of the general analysis of such equations are well emphasized by 
Turing in §11 “Restatement and biological interpretation of the results” (p. 66), 
with an incredible sense of anticipation; it is essential to take into account 

1 . the role of fluctuations, which play a critical role when the system becomes 
unstable; 

2 . the role of slow changes of reaction rates and diffusibility coefficients because 
“such changes are supposed ultimately to bring the system out of the stable 
state”. 


Turing considered therefore that the systems he analyzed belong to the class of 
what are called today slow-fast dynamical systems and focused on the breaking of 
spatial homogeneity near instability. As he said 

“the phenomena when the system is just unstable were the particular 
subject of the inquiry.” 

He underlined the fact that the “linearity assumption” near the equilibrium, i.e. 
the fact that the dynamics is qualitatively equivalent to its linear part, is “a serious 
one” and made what is called today an hypothesis of adiabaticity: as the system is 
a slow-fast one, the slow variation of parameters is slow w.r.t. the fast time used to 
reach equilibrium and, therefore, one can suppose that the system is always in its 
equilibrium state until he reaches a bifurcation destabilizing it. 

In terms of the variables x{r,t) and y{r,t), the linearized react ion-diffusion 
equations are: 


{ x{r, f) = ax{r, f) + by{r, t)-\- 

jj, {x{r — l,t) — 2x{r, t) -|- x{r + 1, t)) 
y{r, t) = cx{r, f) -|- dy{r, t)-|- 

z/ {y{r - 1, t) - 2y{r, t) + y{r + 1, f)) 

In the continuous limit on a circle of radius 1, they are 


l,...,iV) 


x{0,t)\ _ (x{e,t)\ f (x" {6,t) 

y{e,t)) ^^\y{e,t))^\Qv') \y"{e,t) 

where x" and y" are spatial second derivatives (Laplacian term). 

A pedagogical interest of the ring model is that the space is the circle that 
the eigenfunctions of the Laplacian are the trigonometric functions, and that the 
harmonic analysis is therefore nothing else than Fourier analysis. Let f (s, t) and 
f] (s, t) be the Fourier tranforms of x (r, t) and y (r, t)\ 


\v{s,t) 


^Er=f exp(-^)a;(r,t) 
exp {-^)y{r,t) 


(r = l,...,iV) 
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Then x (r, t) and y (r, t) are retrieved through the inverse Fourier transform: 


r a; (r, t) = Yfs=i exp (^) ^ (s, t) 

\ y (f t) = Tfs=i exp (^) r] (s, t) ’ 

By dehnition, the ^ (s, t) and y (s, t) are complex numbers. But as far as x (r, t) and 
y (r, t) are real, we must have ^ (s, t) = ^ {N — s, t) and y (s, t) = y {N — s,t). 

The main interest of using harmonic analysis, is that, in the Fourier domain, the 
system of equations becomes diagonal because the functions are expanded over a 
basis of eigenfunctions of the Laplacian operator. Turing based his computations on 
this separation of variables in the Fourier domain. Due to the dehnition of f (s, t) 

and the expression of x{r,t), the temporal derivatives f {s,t) are 



^{s,t) = 


If one writes rs 
eigenfunctions 


1 

N 


r=N 

X^exp 

r=l 


27iirs\ 

~ir) 


[ax{r, f) + hyij-, t) + y {x{r — 1, t) — 2x{r, f) + x{r- + 1, t))] 

= (r + 1) s — s and uses the orthogonality relations between the 


s=N 


J^expf 

S =1 ^ 


s=N 


277irs\ 

1 = 0 if r = 1,. 

N ) 

27iirs\ 

1 = iV if r = iV 

N J 


then, one gets the equations 


f {s,t) = af (s,t) + by {s,t) + 

/ / 27ris\ f 2Tris\ 

and analog formulae for the y (s, t). One then takes the real and imaginary parts of 
the equations and uses the formulae 


exp 


/ 27ris\ 

\w) 


sm 


2713 


N 


+ sin 


cos 


27rs 

'W 


2 + cos 


2713 

~w 

2713 

~w 


'27rs\ f 2713 

■■ cos I I + I sm 


N 


= 0 


= 2 


cos 


271S 


A ■ 2 

-4 sm^ — 


N 

TIS 

\N 


N 


1 1 = 2 ( cos^ ( ^ j 


since cos 


^l+sm^ 


• 2 

[nJ 

^) = i 

n) 


^ -1 
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to get the equations 


= {a - ifisin^ (f)) ^ {s,t) + bf] {s,t) 

7 ]is,t) = c^{s,t) + {d-Ausin^ (f))r/(s,t) 

In the continuous model, the Fourier transforms of functions on are Fonrier 
series whose components are indexed by /c G Z and one gets 

^ {k, t) = {a- ^ {k, t) + br] {k, t) 

7] {k, t) = {k, t) + {d — T] {k, t) 

with k"^ corresponding to ^ sin^ (^)- Tnring denotes by U this variable. 

5.2 The origin of instability 

The fnndamental new phenomenon introdnced by diffnsion is that the spectral anal¬ 
ysis of the linearized system now depends npon diffusion which, by changing the 
characteristic eqnation, can tranform eigenvalnes with Re (A) < 0 into eigenvalnes 
with Re (A) > 0. It is the origin of diffnsion-driven instabilities. Indeed, the Jacobian 
is now 

/a-d/isin^ (f) b 

V c d-4z/sin2(f) 

and the characteristic eqnation - also called a dispersion relation - is therefore 

— a + Ap, sin^ (^) ) {p ~ dA- Av sin^ (^) ) ^ 

Tnring denotes by ps and with Re (p*) > Re (p^) the two eigenvalnes. If 
Pa 7 ^ p’a then the solntions of the system in the Fonrier domain are of the form 

^ (s, t) = Age^”^ + Baddo^ 

7 ] (s, t) = 

If Pa and p'a are real, then Ajsi-a = Rs, etc. If p* and p'g are conjugate complex 
numbers, then Rtv-s = Ag, etc. It is straightforward to verify that the coefficients 
satisfy the relations: 

( Aa{pa- a + 4/isin^ (^)) = bCg 
1 Bg (p'a -a + Apsin^ (^)) = bDg 

Now if Max (Re {pa)) > 0, some diverging Fonrier modes will become dominant 
and push the system out of equilibrium. Such a possibility can happen only nnder 
precise conditions relating the parameters a, b, c, d of the internal chemical eqnilib- 
rinm to the parameters p, v of the external spatial diffnsion. Tnring explains very 
well that generically only a single Fonrier mode (with its conjngate) can become 
dominant. Indeed, if it was not the case. 
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“the quantities a, 6 , c, d, jj., u will be restricted to satisfy some special 
condition, which they would be unlikely to satisfy by chance.” (p. 50) 

Let So be the index yielding Max (Re (p^)) and suppose Re(psg) > 0. If the two 
eigenvalues Psg and are real, the pair (Pso^Pn-so) will induce divergences since 

sin^ j and if they are complex conjugate the two pairs {pso,PN-so) 

and (Psg,p)v-so) both induce divergences. 

6 A toy model 

Turing presents a simple numerical example p. 52. The parameters are 


a = I-2, 5 = 2.5, c=-1.25, d = J + 1.5 

2 p' V \27ipJ \7ip 


sm 


7rs\ 

ivJ 


The characteristic equation is therefore 


(^p — a + 4p,sin^ (^)) {p ~ 4z/sin^ (^)) ^ 

{p-I+ 2 + U) (^p-I-1.5 + + (2.5) (1.25) = 0 

(P-/M(l + 5f/)(p-/) + t(f/-l)'=o 

We observe that p = I for U = \. Let Sc be the corresponding value of s. If the radius 
p of the ring is such that there exists an integer sq satisfying U = sin^ (q^) ~ 

then there will exist stationary waves with sq lobes. Otherwise, it will be the Sq 
nearest to Sc which will dominate. 

The hguredisplays for / = 0 the graph T of the hyperbola p^+(| + |f/)p + 

I (f/ — 1)^ = 0 in the (17,p) plane for [/ e [0, 1 . 2 ] and p G [—0.4,0], and 1 = 0. 

The points of T are evident. If p is considered as a parameter, 

U = ^ (^1 - 3p ± y/p 2 - lOpj 

and if U is considered as a parameter, 

p = ^ (^-1 - 35/ ± + 1417 - 1 ) . 
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- 0.3 


bi_1_1_1_I_1_1_i_I_1_1_1_1_i_1_1_I_1_1_i_1_i_1_1_u 

0.0 0.2 0.4 0.6 0.8 1.0 1.2 

Figure 3: The graph F of + (| + |?7) p + \ {U — Q ioi U G [0,1.2] and 

p e [-0.4,0]. 

The solutions of + lAU — 1 = 0 are U = —7 ± 5y/2 but, as f/ is a real square 
sin^ (ff)) admissible value is Uc = —7 + 5-\/2 ~ 0.071 and for Uc the 

two values of p are equal to — | (1 + 371) ~ 0.30325. For U > Uc, the two p roots are 
real and for 0 < 7/ < Uc, they have an imaginary part Im (p) = ±-\/77^ + 1477 — 1 
while the real part move on the segment Re (p) = —| (1 + 377) from the point 
(77 = 0,p = —I) to the point (77 ~ 0.071,p ~ 0.30325). In what concerns p, as 77 
is real, we must have -^/p^ — lOp real, that is p^ — lOp > 0 i.e. p G {—oo, 0] or 
p G [10, +oo). 

The figure [^reproduces Turing’s figure 1 which displays Re (p) and — |Im (p)| as 
functions of U for J = 0. 

7 Conditions for instability and the critical point 

In the §9 “Further considerations on the mathematics of the ring”, Turing analyzes 
further the conditions under which a diffusion-driven instability can occur. We will 
present and complete his computations using the continuous model which is easier 
to understand. 

As reaction-diffusion equations are linear, and since every function on is a 
linear superposition of harmonics we look at solutions of the form 

Q = 
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Figure 4: Turing’s figure 1 which displays Re {p) and — |Im (p)| as functions of U for 
1 = 0. The full line and the dotted line represent respectively Re (p*) and Re (p'^), 
while the broken line represents — |Ini(p)|. Tnring has indicated with black thick 
points the integer valnes from s = 0 (left) to s = 5 (right). 

which implies immediately 



The characteristic eqnation is therefore 

Oet (J. - k^D - A/) = Det \ = 0 

and the dispersion relations are 

(a — a + pk'^) (a — d + = be 

or, if we write this eqnation A^ — S' (/c^) A + P (/c^) = 0 with S {k"^) the snm of its 
two roots and P {k^) their prodnet, 

A^ - (Tr (Jo) - k"^ Tr (D)) A + {pu'k^ - {ua + p'd) k'^ + Det (Jo)) = 0 

We want Tr (Jo) = a + d < 0 and Det (Jo) = ad — 6c > 0 to ensnre the stability 
of the internal chemical equilibrinm. Bnt we want also one A with Re (A) > 0 
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to ensure a diffusion-driven instability. As Tr ( Jq) < 0 and Tr (D) = jj,' + v' > 
0, this implies S (/c^) = Tr ( Jq) — Tr {D) <0. If P (/c^) happened to be > 0, 
we would have two roots with Re (A) < 0, so we must have P (/c^) < 0. But as 
Det (Jo) > 0 by hypothesis and of course > 0, we need in fact z/'a -|- jj,'d > 

-|- Det (Jo) /k"^ > 0. This is a first condition. As Tr (Jo) = a -f- J < 0 and 
/r', P > 0 , we need n' ^ P since, if n' and v' would be equal, we would have 
z/'a -|- n'd = jj! {a + d) < 0 . 

It is essential to strongly emphasize here the fact that the diffusibility of the two 
morphogens X and Y must be sufficiently different in order that an instability can 
occur. It is the key of Turing’s discovery. 

In the example of §10, /i' = q;|, z/' = q;| (a > 0), a = —2, J = | and the 
condition z/'a -|- fi'd = a (—2 x|-|-|x^) = |> 0 is therefore satisfied. 

There is a second condition. P (/c^) = /i'z/'/c^ — (z/'a -|- /r'J) k'^ + Det (Jo) is a second 
degree polynomial in k‘^ and we want it to become < 0 for values of k"^ which must 
necessarilly be positive since is a real square. Let 6 = ^. The graph of P (A;^) is 
a parabola 11 starting at Det (Jo) > 0 for A; = 0. If 5 < 5c for a critical value to be 
computed, 11 is over the A;^-axis and the condition P (A;^) < 0 cannot be satisfied. 
But if 5 > Sc, n intersects the A;^-axis at two points kl and k"^ > kf and inside the 
interval [kf,kl] the condition P {k"^) < 0 is satisfied: there exists an eigenvalue A 
with Re (A) > 0. 

The computation of Sc is rather tedious. Let u = k‘^. The polynomial P {u) and 
its first and second derivatives are 


P (m) = — (z/'a -|- fi'd) u + ad — be 

P' (u) = — (z/'a -|- /u'd) + 2n'ffu 
P" {u)=2fiff' 


So the minimum of 11 is given by P' (m) = — (z/'a -|- /a'd) + 2ij,'n'u = 0 and it is a 
minimum since P" {u) = 2/i'z/' > 0. Its value is 


and we have 


uo = kl 


z/'a -|- fj,'d 
2/i'z/' 


P{kl)=ad-bc-tiY±lPl 

^ ^ 4 /i'z/' 

For P (A;^) to become < 0, we must have therefore 

, , 1 (z/'a -|- /i'J)^ 1 (a -|- SdY 

0<ad-bc< ^ 

4 /i z/ 4 S 

and Sc is given by the equation 
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ad — be 


1 (a + 5cd)^ 

4 5 ;^ 

In the example, when 7 = 0, we have ad — be = u'a + jj!d = /iV' = and 
we verify that | = | x (|)^ x 8 . So Turing’s system is at its critical point for 7 = 0 . 

Let us now investigate more precisely the second degree equation giving Sc- It 
can be written 


d‘^5l + 2 (26c — ad) 5c + a"^ = 0 

and its two solutions are 

ad — 2bc ± 2 \/—bc {ad — be) 

But as one root, and hence both roots, must be real, we need —be {ad — be) > 0 and, 
as ad — 6 c > 0 by hypothesis, we must have 6 c < 0 . 

In the example, for 7 = 0, we have effectively —x 2 < 0 and 



As 6 = 6c we are indeed at the critical point. Then 

and the parabola 11 is tangent at the /c^-axis and at the point of tangency the 
eigenvalues are A+ = 0 and A_ = Tr ( Jq) — k'^ Tr {D) = — — | < 0. So, at the 

crossing of the critical point, the eigenvalue A+ becomes > 0. 

8 The bifurcation 

After having analyzed the conditions of a diffusion-driven instability at a critical 
point, Turing analyzed further the behavior of the system in the neighbourhood of 
the critical point. To this end, he varied the small parameter 7 around its critical 
value 7 = 0 . Today, computations are very easy, but at his time they were difficult 
and he must use the computer he had himself constructed. We will hrst do them 
for the continuous model and then return to Turing’s own discrete model. 

8.1 Continuous model 

The chemical internal equilibrium is now given by the concentrations of morphogens 

Xo = ^ (-25 + v /210 + 7 X 557 ) , Xo = 1 
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Figure 5: The graph of P (/c^) for e = 0.1. Inside the interval [0.367, 2.862] of we 
have P < 0. 


We linearize the system in the neighbourhood of (Xq, To) and compute hrst order 
expansions in the small parameter £ = -^11^7 = ^^^7- We get 


a = —2 — e,b 


25 

16 


25 

—e, c = 2 + e,d 

7 



and the conditions for instability: Tr ( Jq) = a + d = —1 + must be < 0, which 
implies £ < ~ 0.194; Pa + n'd = | + must be > 0, which implies e > — 

Det (Jo) = ad — be = | + must be > 0, which implies e > —2; ad must be < 0, 

which implies e > ad — be must be < | , which implies £ > 0. 

The equation yielding the eigenvalues is now 


- (Tr (Jo) - P Tr (D)) X + P {k^) = 0 

Figure 1^ displays the graph of P (iP) for e = 0.1. The roots of P (IP) = 0 are 
0.367 and 2.862 and inside their interval we have P (fc^) < 0. 

The characteristic equation is 


- (Tr (Jo) - Tr {D)) A + (/iV'fc^ - {Pa + /I'J) k"^ + Det (Jo)) = 0 

that is 


A2 + 


fl 3fc2 

(2 + x 


18 . , 

1 + 



and its roots A± are approximated by 






0 


1 


(^-7 + 36e ± \/-49 - 553£ + 1296£2 + 490fc2 - 3085^2 + 343^4^ 


Figures and display the graphs of A+ and A (including the irrelevant negative 
fc^-axis). For A+ we see that A+ > 0 for G [0.367,2.862]. Figure zooms on 
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Figure 6: The discriminant A for £ = 0.1. A is negative inside the open interval 
]-1.5146,0.1758[. 



Figure 7: Graph of A+ (including the irrelevant negative fc^-axis). A+ > 0 for 
k‘^ G [0.367,2.862]. There is no graph inside the open interval ] —1.5146, 0.1758[ 
where the discriminant A of the characteristic equation is < 0. 


this interval. There is no graph inside the open interval ]—1.5146, 0.1758[ where the 
discriminant A of the characteristic equation is < 0 (see figure]^. 

8.2 Discrete model 

Let us come back to the discrete ring model composed of A = 20 cells. At the 
critical point 7 = 0, the 20 characteristic equations in the Fourier domain are 

(p + 2 + 2 W (1^)) (p - 1.5 + (I)) + I = 0 

The hgure[l0| shows the table of the 20 pairs {ps,p's) of eigenvalues for s = 0,..., 19. 
We see that if we order the p w.r.t to increasing Re (p) we get ps = pn = —0.00346, 
Pa = Pw = —0.012, ps = pi5 = —0.064, p2 = pis = —0.066. It is therefore the mode 
P3 which can most readily become > 0. 
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Figure 8: Zoom on the interval G [0.367,2.862] of the hgure[^ where A+ > 0. 



Figure 9: Graph of A_ (including the irrelevant negative fc^-axis). A_ is always < 0 
for > 0. There is no graph inside the open interval ]—1.5146, 0.1758[ where the 
discriminant A of the characteristic equation is < 0. 


eq [ s_] := (p + 2 + 2* (Sin[Pi*s/20])''2) * (p - 1.5 + (Sin[Pi*s/20])''2) + 25/8 == 0; 

Table[Solve[eq[s], p], {s, 0, 19}] 

{{{p -> -0.25 - 0.251}, {p -> -0.25 + 0.251}}, 

{{p -> -0.286708 - 0.1397311}, {p -> -0.286708 + 0.1397311}}, 

{{p -> -0.720177}, {p -> -0.0662972}}, 

{{p -> -1.11487}, {p -> -0.00345613}}, {{p -> -1.52451}, {p -> -0.0119627}}, 

{{p -> -1.93541}, {p -> -0.0645857}}, {{p -> -2.32263}, {p -> -0.140898}}, 

{{p -> -2.65919}, {p -> -0.222488}}, {{p -> -2.92013}, {p -> -0.293399}}, 

{{p -> -3.08537}, {p -> -0.341218}}, {{p -> -3.14194}, {p -> -0.358059}}, 

{{p -> -3.08537}, {p -> -0.341218}}, {{p -> -2.92013}, {p -> -0.293399}}, 

{{p -> -2.65919}, {p -> -0.222488}}, {{p -> -2.32263}, {p -> -0.140898}}, 

{{p -> -1.93541}, {p -> -0.0645857}}, {{p -> -1.52451}, {p -> -0.0119627}}, 

{{p -> -1.11487}, {p -> -0.00345613}}, {{p -> -0.720177}, {p -> -0.0662972}}, 

{{p -> -0.286708 - 0.139731 1}, {p -> -0.286708 + 0.139731 1}}} 

Figure 10: The table of the 20 pairs {ps.v's) of eigenvalues of the discrete ring model 
for s = 0,..., 19 and 7 = 0. 
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eq2[s_] := 

{p + 2.02336 + 2* (Sin [Pi *s/20])-'2) * (p-1.58344 + (Sin [Pi *s/20])-'2) ==-1.64594*2.02 336 
iist2 = Table[Solve[eq2[s], p], {s, 0, 19}] 

{{{p^-0.21996-0.2794241}, (p^-O.21996 + 0.2794241}}, 

{{p ^ - 0.256668 - 0.183836 1}, {p ^ - 0.256668 + 0.183836 1}}, 

{ {p ^ - 0.673699}, (p ^ -0.0526953}}, 

{{p ^ - 1.0807}, {p ^ 0.0224553}}, {{p ^ - 1.49637}, {p ^ 0.0199735}}, 

{ {p ^ -1.9113}, {p ^ - 0.0286193}}, {{p ^ - 2.30143}, (p ^ - 0.102014}}, 

{ {p ^ -2.64011}, (p ^ - 0.181492}}, {{p ^ - 2.90249}, (p ^ - 0.25096}}, 

{{p ^ - 3.06857}, {p ^ - 0.297935}}, {{p ^ - 3.12542}, {p ^ - 0.314498}}, 

{{p ^ - 3.06857}, {p ^ - 0.297935}}, {{p ^ - 2.90249}, (p ^ - 0.25096}}, 

{{p-» -2.64011}, {p ^ - 0.181492}}, {{p-» - 2.30143}, {p ^ - 0.102014}}, 

{ {p ^ -1.9113}, {p ^ - 0.0286193}}, {(p ^ - 1.49637}, (p ^ 0.0199735}}, 

{{p ^ - 1.0807}, {p ^ 0.0224553}}, {{p ^ - 0.673699}, {p ^ - 0.0526953}}, 

{{p ^ - 0.256668 - 0.183836 i}, {p ^ - 0.256668 + 0.183836 i}}} 


Figure 11: The table of the 20 pairs {ps,p's) of eigenvalues of the discrete ring model 
for s = 0,..., 19 and 1 = 


Let us now vary the small slow parameter 7. Turing varied 7 almost adiabatically 
from — f (stability) to (instability) at speed 7 = 2“^ = This corresponds 
to variations e : —0.094 —)■ 0.0235 for e and f : 0 —?■ 40 for the discrete time t. 
For 7 = —f, the equilibrium is {Xq = 0.78,Yq = I), and Oq = —1.9, Bq = —1.218, 
Co = 1.9, do = 1.156, and all Re (ps) < 0: the system is stable. On the contrary, for 
7 = —^, the equilibrium is {Xi = 1.053, Yi = 1), and Oi = —2.023, hi = —1.646, 
Cl = 2.023, di = 1.583, and ps = pi7 = 0.0224 > 0, p4 = pio = 0.012 > 0: the system 
is unstable. The hgure 11 shows the table of the 20 pairs {ps,Ps) of eigenvalues for 
s = 0,..., 19 for 7 = ^. 

In figure 12 we show the 20 graphs Re {ps) as functions of t for 7 varying from 
— f to The time t varies from f = 0 to f = 40. We see the Re {ps) which become 
> 0 for s = 3,17,4,16. At f = 40, ps = pn = 0.0224 and p4 = pie = 0.012. For 
s = 0,1,19, Re {ps) presents an angular point because ps is a complex number with 
Im(ps) 7^ 0. It is the same phenomenon as in the toy model of figure]^ Figure 13 
shows the graphs of Re (p) and Im (p) in such a case. 

In his paper, Turing computes (with his recently constructed Manchester com¬ 
puter) the table of the evolution of the ring (see hgure [l^ and shows how Fourier 
modes become dominant after the bifurcation induced by the diffusion-driven insta¬ 
bility (see hgure 15). In the initial state, all cells are, up to small huctuations, in 
the equilibrium state (Xq = 1,1o = !)• After the bifurcation, a stationary oscilla¬ 
tory wave pattern with 3 lobes develops. The divergences induced by the instability 
are tamed by two factors: (i) the concentration X cannot become < 0 and when 
X (r, t) vanishes, the process stops locally^ (ii) saturation non-linear ehects allow a 
new equilibrium to occur. These results constitute a great achievement. 


9 Further aspects of Turing’s paper 

In his paper, Turing evoked many other problems. In §4 he gave simple examples 
for explaining the idea of “breakdown of symmetry and homogeneity” in pattern 
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Figure 12: The 20 graphs Re {ps) as functions of t for 7 varying from —| to The 
time t varies from t = 0 to t = 40. The Re {ps) become > 0 for s = 3,17,4,16. At 
t = 40, Ps = pi 7 = 0.0224 and p4 = piQ = 0.012. For s = 0,1,19, Re (ps) presents 
an angular point because ps is a complex number with Im (ps) 7^ 0. It is the same 
phenomenon as in the toy model of hgure 
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Figure 13: When Im (p) 7^ 0 (graph up), Re (p) (graph down) presents an angular 
point. 


Table 1. Some stationary-wave patterns 




first 

specimen 


cell 

incipient pattern 

final pattern 

number 

A' 

Y 

X 

Y 

0 

1130 

0-929 

0-741 

1-463 

1 

1123 

0-940 

0-761 

1-469 

2 

1164 

0-885 

0-954 

1-255 

3 

l-21o 

0-810 

1-711 

0-0(K1 

4 

1-249 

0-753 

1-707 

0-000 

5 

1-158 

0-873 

0-875 

1-385 

6 

1-074 

1-003 

0-700 

1-622 

7 

1-078 

1-000 

0-699 

1-615 

8 

1-148 

0-896 

0-885 

1-382 

9 

1-231 

0-775 

1-704 

0-000 

10 

1-204 

0-820 

1-708 

0-000 

11 

1-149 

0-907 

0-944 

1-273 

12 

1-156 

0-886 

0-766 

1-461 

13 

1-170 

0-854 

0-744 

1-442 

14 

1-131 

0-904 

0-756 

1-478 

15 

1-090 

0-976 

0-935 

1-308 

16 

1-109 

0-957 

1-711 

0-000 

17 

1-201 

0-820 

1-706 

0-000 

18 

1-306 

0-675 

0-927 

1-309 

19 

1-217 

0-811 

0-746 

1-487 


second 

‘slow 



specimen: 

incipient 

Y 

cooking’: 

incipient 

Y 

four-lobed equilibrium 

X 

Y 

0-834 

1-057 

1-747 

0-000 

0-833 

0-903 

1-686 

0-000 

0-766 

0-813 

1-445 

2-500 

0-836 

0-882 

0-445 

2-500 

0-9:10 

1-088 

1-685 

0-000 

0-898 

1-222 

1-747 

0-000 

0-770 

1-173 

1-685 

0-000 

0-740 

0-956 

0-445 

2-500 

0-846 

0-775 

0-445 

2-500 

0-937 

0-775 

1-685 

0-000 

0-986 

0-969 

1-747 

0-000 

1-019 

1-170 

1-685 

0-000 

0-899 

1-203 

0-445 

2-500 

0-431 

1-048 

0-445 

2-500 

0-485 

0-868 

1-685 

0-000 

0-919 

0-813 

1-747 

0-000 

1-035 

0-910 

1-685 

0-000 

1-003 

1-050 

0-445 

2-500 

0-899 

1-175 

0-445 

2-500 

0-820 

1-181 

1-685 

0-000 


Figure 14: Turing’s computation of the evolution of the ring. 
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Figure 15: The evolution of the ring after the bifurcation induced by a diffusion- 
driven instability. The hatched graph represents the X concentration at the initial 
state: all cells are, up to small fluctuations, in the equilibrium state (Xq = 1, Tq = !)• 
The other graph represents a stationary oscillatory wave pattern with 3 lobes. The 
divergences induced by the stability are tamed by two factors: (i) the concentration 
X cannot become < 0, (ii) non-linear effects. 
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Figure 16: Turing’s figure 2 on “dappled colour patterns” in two-dimensionsional 
tissues. 

formation. He explained also how exponential divergences are bounded by non- 
linearities which allow new equilibria to emerge and he emphasized 

“the effect of considering non-linear reaction rate functions when far from 
homogeneity.” (p. 58) 

In §8 he listed some “types of asymptotic behaviours in the ring after a lapse of 
time”: 

1. “stationary cases” where the asymptotic regime is dominated by a pair of real 
eigenvalues (pso^Pso)'^ 

2. “oscillatory cases” where the asymptotic regime is dominated by a pair of 
complex conjugated eigenvalues {pso,p'so) (travelling waves); 

3. “limit cases”. 

He drew also phase diagrams in the parameter space which classify these different 
regimes and explained the role of fluctuations in the bifurcation process. 

Another extremely important anticipation is that in two-dimensionsional tissues 
diffusion-driven instabilities can explain “dappled colour patterns” as observed on 
sea shells or leopard’s coats. Figure [16] reproduces Turing’s hgure 2. 

Moreover, Turing anticipated the fact that his general model was able to induce 
oscillating patterns when the chemical internal dynamics of each cell bifurcates to¬ 
wards a limit cycle by Hopf bifurcation. When such limit cycles propagate spatially, 
many complex phenomena can emerge. Turing envisaged applications to organisms 
such as plants (flowers, leaves) or Hydra. His predictions have been widely conhrmed 
later. 

In the fascinating §12 “Chemical waves on spheres. Gastrulation”, Turing gen¬ 
eralizes his one-dimensional ring model to a two-dimensional sphere model whose 
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geometry is more complex, the harmonic analysis on the sphere resting on the eigen¬ 
functions of the spherical Laplacian, namely the spherical harmonics. His striking 
idea was to apply the model to gastrulation in embryology, which is the step at which 
the spherical symmetry of the blastula is broken. He developed this idea further in 
his unpublished paper on phyllotaxis. 

Finally, in the last §13 “Non-linear theory. Use of digital computer”, Turing came 
back to the assumption that linearization is a good approximation and explained 
that it is the case only in the neighbourhood of the hrst bifurcation. It is 

“an assumption which is justihable in the case of a system just beginning 
to leave a homogeneous condition.” (p. 72) 

For Turing it was risky to try to go beyond: 

“One cannot hope to have any very embracing theory of such processes, 
beyond the statement of the equations.” (p. 73) 

Hence the fundamental interest of the just constructed digital computers enabling 
numerical simulations avoiding the too drastic simplihcations imposed by the search 
of explicit theoretical solutions to the equations. 

10 Conclusion: after Turing 

To conclude this presentation of Turing’s 1952 paper, let us look briefly at the works 
on reaction-diffusion equations after Turing. I have already tackle this theme in my 
talk [19] at the TAPS 2001 Conference on Complexity and Emergence. 

10.1 General reaction-diffusion models 

Among the many specialists of the domain, we would cite Hans Meinhardt and 
Alfred Gierer who, since 1972 |S|, have considerably increased our knowledge on 
reaction-diffusion models. They have shown that, for an activator/inhibitor pair of 
morphogens, instabilities mainly result from the competition between a short range 
slow activation and a long range fast inhibition, the inhibitor diffusing faster than 
the excitator. This conhrms Turing’s remark on the role of the difference between 
the two diffusibility coefficients /i and u. 

A general Meinhardt-Gierer model has the form 

( X = p— — ax + ax + ZiAx „ 

{ . ^ n A a<P, 

[y = px^ - (3y + ay + nAy 

The activator morphogen x is self-catalizing term in x) and its production is 
inhibited by the inhibitor morphogen y term in x). Moreover, x catalyzes its 
inhibitor term in y). The linear terms ax and /3y {a < (3) are degradation terms, 
the constant ay enables a stable homogeneous state and the constant a^ allows to 
trigger the process, p and n are the diffusibility coefficients with p (slow) -C ly 
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Figure 17: The coupling between internal dynamics and external diffusion can induce 
very complex patterns. Two examples due to De Kepper: a system of stripes (with 
defects) and a honeycomb pattern. (From De Kepper et al. jSj). 

(fast). A local fluctuation of the activator x induces a local peak of x which diffuses 
slowly. But it amplihes also the inhibitor concentration y, and since y diffuses faster 
than X it will inhibit the production of x at some distance (what is called a “lateral 
inhibition”). 

The coupling between the internal dynamics and external diffusion can induce 
very complex patterns. Figure [l7| shows two examples due to De Kepper [3], a system 
of stripes with defects and a honeycomb pattern. 

One of the best known achievements of Hans Meinhardt is his modelling of sea 
shells. Since the growth of a shell results from an accretion of calcified matter 
along its boundary, it can be represented by a two-dimensional diagram B x T = 
boundary X time. The geometry is therefore in fact one-dimensional. The diffusion 
of the activator from a local peak of concentration induces a triangle where the 
pigmentation controlled by x is high, but the faster diffusion of the inhibitor stops 
it after a while. Hence a cascade of triangles. Figure [TS] shows the celebrated model 
of Conus marmoreous. 

10.2 Prom Alan Turing to Rene Thom 

As we have seen, in Turing’s paper morphogenetic processes spatially unfold diffusion- 
driven instabilities. In the late 1960s, Rene Thom I2D1, IZH proposed a more general 
model based on the general concept of bifurcation. The similarities and dissimilari¬ 
ties between Turing’s and Thom’s models are fascinating. The key idea is the same: 
internal chimical dynamics (reactions) are coupled with external spatial dynamics, 
the latter destabilize the former and morphologies spatially unfold the instabilities. 
As we have seen, in Turing coupling is given by the diffusion of the reacting mor- 
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Figure 18: Meinhardt’s model for the sea shell Conus marmoreous. In front a true 
shell. In the background its reaction-diffusion model. (From Meinhardt [TU|h 


phogens. In Thom, coupling is more generally a spatial control of internal dynamics 
and the morphogenetic discontinuities breaking the homogeneity of the substrate 
are induced by bifurcations. 

10.3 Beyond Turing 

After Turing, many authors, e.g. James Murray na, na, introduced bifurcations 
in reaction-diffusion equations. Let u be the vector (x, y) and consider a differential 
equation u = f {u, r) where r is a spatial control. When r varies and crosses a critical 
value Tc, the initial stable equilibrium state Uq of the system can collapse with an 
unstable equilibrium and disappear. The system is therefore projected to another 
equilibrium through this saddle-node bifurcation. Another most used bifurcation 
is the Hopf bifurcation. When r varies and crosses a critical value rc, the initial 
stable (i.e. attracting) equilibrium state uq becomes a repellor and generates a 
small attracting closed orbit (i.e. a limit cycle). 

Consider for instance the following system analyzed by Robin Engelhard! |1]: 

j X = —xy"^ + ay — {1 + h) X + 6Ax 
1 1 / = xy‘^ — {l-\-a)y-\-x-\-F-\- 6 Ay 

The chemical internal equilibria without diffusion {6 = 0) are solutions of the equa¬ 
tions (if 1 /^ J- 1 J- 6 7 ^ 0): 
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that is 


X = 


ay 


+ 1 + b 


y 


y‘^ + 1 + b 


{l + a)y + 


ay 


y'^ + 1 + b 


+ F = 0 


y^ — Fy"^ + {1 + b + ab) y — F {1 + b) = 0 

which is a cubic equation with parameters a, 6, F. 

At the points where ?/^ + 1 + 6 = 0, we have 

j X = ay + 6 Ax 

I ^ = —bx — {1 + a)y + F + 6 Ay 

and this can be an equilibrium point for 5 = 0 only ii ay = 0 and —bx — y + F = 0. 
If a 7 ^ 0, this implies the condition b = —1, and the equilibrium is ?/ = 0, a; = —F. 


shows some examples of patterns solution of this system of equations. 

There is a wealth of material on these topics. The reader could look e.g. at 
Harrison [U], Lee et al. 0 . 0 . Maini [U], Oyang-Swinney [IB], or Pearson m 

10.4 Experimental results 

The validity of Turing’s models for embryology are still under discussion. But in 
what concerns chimical and biological patterns their validity is without doubt. We 
have seen Meinhardt’s examples. For chemical systems exact verihcations go back 
to 1990 and the works of the Bordeaux group of Patrick De Kepper (Castets, Du- 
los, Boissonade) on iodate-ferrocyanide-sulhte or clorite-iodide-malonic acid-starch 
reactions in gel reactors. 

It is a full universe of morphological phenomena and mathematical models that 
Turing opened in 1952 with a remarkable foresightedness. 
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Figure 19: Some examples of patternized solutions of Engelhardt’s system of equa¬ 
tions. (From |1]). 
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